//==============================================================================
//==============================================================================
/// \file stellarpop_data.hpp
/// \brief Data to extract fate and lifetime of high mass stars
///

// #ifndef PLANCKINTEGRAL_HPP_ // NOLINT
// #define PLANCKINTEGRAL_HPP_

#include <cassert>
#include <cmath>

#include "AMReX_Array.H"
#include "AMReX_BLassert.H"
#include "AMReX_Extension.H"
#include "AMReX_GpuQualifiers.H"
#include "AMReX_REAL.H"
#include "fundamental_constants.H"
#include "math/interpolate.hpp"

using Real = amrex::Real;
static constexpr int FATE_ARR_SIZE = 201;
static constexpr int AGE_ARR_SIZE = 197;
static constexpr int YR_TO_SEC = 3.15576e7; // seconds in a year

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto interpolate_fate(Real mass_star) -> int
{
	const amrex::GpuArray<Real, FATE_ARR_SIZE> interp_mass_star = {
	    // mass of stars in Msun
	    0.0,  9.,	9.25, 9.5,  9.75, 10.,	10.25, 10.5, 10.75, 11.,  11.25, 11.5, 11.75, 12.,  12.25, 12.5, 12.75, 13.,  13.1, 13.2, 13.3, 13.4, 13.5,
	    13.6, 13.7, 13.8, 13.9, 14.,  14.1, 14.2,  14.3, 14.4,  14.5, 14.6,	 14.7, 14.8,  14.9, 15.,   15.1, 15.2,	15.3, 15.4, 15.5, 15.6, 15.7, 15.8,
	    15.9, 16.,	16.1, 16.2, 16.3, 16.4, 16.5,  16.6, 16.7,  16.8, 16.9,	 17.,  17.1,  17.2, 17.3,  17.4, 17.5,	17.6, 17.7, 17.8, 17.9, 18.,  18.1,
	    18.2, 18.3, 18.4, 18.5, 18.6, 18.7, 18.8,  18.9, 19.,   19.1, 19.2,	 19.3, 19.4,  19.5, 19.6,  19.7, 19.8,	19.9, 20.,  20.1, 20.2, 20.3, 20.4,
	    20.5, 20.6, 20.7, 20.8, 20.9, 21.,	21.1,  21.2, 21.3,  21.4, 21.5,	 21.6, 21.7,  21.8, 21.9,  22.,	 22.1,	22.2, 22.3, 22.4, 22.5, 22.6, 22.7,
	    22.8, 22.9, 23.,  23.1, 23.2, 23.3, 23.4,  23.5, 23.6,  23.7, 23.8,	 23.9, 24.,   24.1, 24.2,  24.3, 24.4,	24.5, 24.6, 24.7, 24.8, 24.9, 25.,
	    25.1, 25.2, 25.3, 25.4, 25.5, 25.6, 25.7,  25.8, 25.9,  26.,  26.1,	 26.2, 26.3,  26.4, 26.5,  26.6, 26.7,	26.8, 26.9, 27.,  27.1, 27.2, 27.3,
	    27.4, 27.5, 27.6, 27.7, 27.8, 27.9, 28.,   28.1, 28.2,  28.3, 28.4,	 28.5, 28.6,  28.7, 28.8,  28.9, 29.,	29.1, 29.2, 29.3, 29.4, 29.5, 29.6,
	    29.7, 29.8, 29.9, 30.,  31.,  32.,	33.,   35.,  40.,   45.,  50.,	 55.,  60.,   70.,  80.,   100., 120.};

	const amrex::GpuArray<Real, FATE_ARR_SIZE> interp_star_fate = {
	    // if fate==1, star will go SN, inject energy!
	    // if fate==0, star will not go SN, no energy injection
	    0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
	    1., 1., 1., 0., 0., 1., 0., 0., 0., 0., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 0., 0., 1., 0.,
	    1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 0., 0., 0., 1., 1., 0., 0., 1., 1., 1., 1., 1., 0., 0., 1., 0., 1., 1., 1., 0., 0.,
	    1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
	    0., 0., 0., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 0., 0., 0., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
	    0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 1.};

	// Interpolate to find the accurate fates from array
	auto const &x_arr = interp_mass_star;
	auto const &y_arr = interp_star_fate;
	const double mass_in_Msun = mass_star / C::M_solar;
	AMREX_ASSERT(mass_in_Msun >= 0.0);

	// Interpolate to find the fate of all masses, using clamp policy to return first/last element for out-of-bounds
	amrex::Real fate_interp = interpolate_value<BoundaryPolicy::Clamp>(mass_in_Msun, x_arr.data(), y_arr.data(), FATE_ARR_SIZE); // NOLINT
	return (fate_interp < 0.5 ? 0 : 1);
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto interpolate_death_time(Real mass_star) -> Real
{
	const amrex::GpuArray<Real, AGE_ARR_SIZE> interp_mass_star = {
	    // mass of stars in Msun
	    0.0,      1.00e-01, 1.50e-01, 2.00e-01, 2.50e-01, 3.00e-01, 3.10e-01, 3.20e-01, 3.30e-01, 3.40e-01, 3.50e-01, 3.60e-01, 3.70e-01, 3.80e-01,
	    3.90e-01, 4.00e-01, 4.50e-01, 5.00e-01, 5.50e-01, 6.00e-01, 6.50e-01, 7.00e-01, 7.50e-01, 8.00e-01, 8.50e-01, 9.00e-01, 9.20e-01, 9.40e-01,
	    9.60e-01, 9.80e-01, 1.00e+00, 1.02e+00, 1.04e+00, 1.06e+00, 1.08e+00, 1.10e+00, 1.12e+00, 1.14e+00, 1.16e+00, 1.18e+00, 1.20e+00, 1.22e+00,
	    1.24e+00, 1.26e+00, 1.28e+00, 1.30e+00, 1.32e+00, 1.34e+00, 1.36e+00, 1.38e+00, 1.40e+00, 1.42e+00, 1.44e+00, 1.46e+00, 1.48e+00, 1.50e+00,
	    1.52e+00, 1.54e+00, 1.56e+00, 1.58e+00, 1.60e+00, 1.62e+00, 1.64e+00, 1.66e+00, 1.68e+00, 1.70e+00, 1.72e+00, 1.74e+00, 1.76e+00, 1.78e+00,
	    1.80e+00, 1.82e+00, 1.84e+00, 1.86e+00, 1.88e+00, 1.90e+00, 1.92e+00, 1.94e+00, 1.96e+00, 1.98e+00, 2.00e+00, 2.02e+00, 2.04e+00, 2.06e+00,
	    2.08e+00, 2.10e+00, 2.12e+00, 2.14e+00, 2.16e+00, 2.18e+00, 2.20e+00, 2.22e+00, 2.24e+00, 2.26e+00, 2.28e+00, 2.30e+00, 2.32e+00, 2.34e+00,
	    2.36e+00, 2.38e+00, 2.40e+00, 2.42e+00, 2.44e+00, 2.46e+00, 2.48e+00, 2.50e+00, 2.52e+00, 2.54e+00, 2.56e+00, 2.58e+00, 2.60e+00, 2.62e+00,
	    2.64e+00, 2.66e+00, 2.68e+00, 2.70e+00, 2.72e+00, 2.74e+00, 2.76e+00, 2.78e+00, 2.80e+00, 3.00e+00, 3.20e+00, 3.40e+00, 3.60e+00, 3.80e+00,
	    4.00e+00, 4.20e+00, 4.40e+00, 4.60e+00, 4.80e+00, 5.00e+00, 5.20e+00, 5.40e+00, 5.60e+00, 5.80e+00, 6.00e+00, 6.20e+00, 6.40e+00, 6.60e+00,
	    6.80e+00, 7.00e+00, 7.20e+00, 7.40e+00, 7.60e+00, 7.80e+00, 8.00e+00, 9.00e+00, 1.00e+01, 1.10e+01, 1.20e+01, 1.30e+01, 1.40e+01, 1.50e+01,
	    1.60e+01, 1.70e+01, 1.80e+01, 1.90e+01, 2.00e+01, 2.20e+01, 2.40e+01, 2.60e+01, 2.80e+01, 3.00e+01, 3.20e+01, 3.40e+01, 3.60e+01, 3.80e+01,
	    4.00e+01, 4.50e+01, 5.00e+01, 5.50e+01, 6.00e+01, 6.50e+01, 7.00e+01, 7.50e+01, 8.00e+01, 8.50e+01, 9.00e+01, 9.50e+01, 1.00e+02, 1.05e+02,
	    1.10e+02, 1.15e+02, 1.20e+02, 1.25e+02, 1.30e+02, 1.35e+02, 1.40e+02, 1.45e+02, 1.50e+02, 1.75e+02, 2.00e+02, 2.25e+02, 2.50e+02, 2.75e+02,
	    3.00e+02};

	const amrex::GpuArray<Real, AGE_ARR_SIZE> interp_death_time{
	    // Age of star upon death in yr
	    // quick hack: assume the age of stars below 0.1 Msun is the same as that of 0.1 Msun. To avoid out-of-bounds interpolation.
	    2.05637091e+12, 2.05637091e+12, 1.23431079e+12, 8.04320575e+11, 5.28141713e+11, 3.23462726e+11, 2.90832810e+11, 2.50741129e+11, 2.19855225e+11,
	    1.91893537e+11, 1.58212735e+11, 1.47856310e+11, 1.39272077e+11, 1.31394380e+11, 1.24030878e+11, 1.17112455e+11, 8.65504379e+10, 6.11287818e+10,
	    4.38991462e+10, 3.40965212e+10, 2.58942710e+10, 2.00010752e+10, 1.71445535e+10, 1.37414072e+10, 1.11618939e+10, 9.18190797e+09, 8.52317286e+09,
	    7.92367923e+09, 7.38317517e+09, 6.88750285e+09, 6.43787571e+09, 6.02635418e+09, 5.65256243e+09, 5.30792038e+09, 4.99571866e+09, 4.70398362e+09,
	    4.43604634e+09, 4.19160550e+09, 3.95903488e+09, 3.74780723e+09, 3.55280230e+09, 3.37061965e+09, 3.20145766e+09, 3.04574464e+09, 2.90039962e+09,
	    2.76673510e+09, 2.64070528e+09, 2.52358989e+09, 2.41321395e+09, 2.31170278e+09, 2.21767859e+09, 2.12845069e+09, 2.04768934e+09, 1.96795341e+09,
	    1.89655450e+09, 1.82697024e+09, 1.76205837e+09, 1.70319680e+09, 1.64588557e+09, 1.59304461e+09, 1.54180672e+09, 1.49204275e+09, 1.44639718e+09,
	    1.40122675e+09, 1.35912000e+09, 1.31910886e+09, 1.28399834e+09, 1.24729434e+09, 1.21422515e+09, 1.18132902e+09, 1.15638835e+09, 1.12748058e+09,
	    1.10297869e+09, 1.08638131e+09, 1.08204826e+09, 1.07007878e+09, 1.04709408e+09, 1.02263616e+09, 9.98411328e+08, 9.74186560e+08, 9.49961792e+08,
	    9.23821184e+08, 9.00427520e+08, 8.73792640e+08, 8.53875648e+08, 8.28805056e+08, 8.08528192e+08, 7.87503552e+08, 7.68840576e+08, 7.48485760e+08,
	    7.30946880e+08, 7.13422336e+08, 6.94482816e+08, 6.80367744e+08, 6.63504000e+08, 6.48566400e+08, 6.33002048e+08, 6.19533120e+08, 6.05064256e+08,
	    5.92603904e+08, 5.77887872e+08, 5.66185856e+08, 5.53863360e+08, 5.43793408e+08, 5.30384992e+08, 5.20402528e+08, 5.07376320e+08, 4.99770336e+08,
	    4.88929568e+08, 4.79426656e+08, 4.69923776e+08, 4.60959264e+08, 4.52413152e+08, 4.42090624e+08, 4.35036448e+08, 4.27982272e+08, 4.19835296e+08,
	    4.11688320e+08, 4.03461568e+08, 3.95234848e+08, 3.88593568e+08, 3.26886272e+08, 2.80325856e+08, 2.45220128e+08, 2.12310240e+08, 1.88788672e+08,
	    1.67970992e+08, 1.51119456e+08, 1.36743856e+08, 1.24152872e+08, 1.13495016e+08, 1.03689080e+08, 9.43302000e+07, 8.74026240e+07, 8.11816960e+07,
	    7.51093840e+07, 6.97079920e+07, 6.62092000e+07, 6.27928920e+07, 5.93765840e+07, 5.59602760e+07, 5.25439680e+07, 4.91276600e+07, 4.67201800e+07,
	    4.39142720e+07, 4.21476160e+07, 4.03531400e+07, 3.32273580e+07, 2.65855740e+07, 2.27899720e+07, 1.99802960e+07, 1.78570140e+07, 1.61082790e+07,
	    1.46735110e+07, 1.35130840e+07, 1.25814310e+07, 1.17604370e+07, 1.10608340e+07, 1.04301560e+07, 9.43675700e+06, 8.75512700e+06, 8.42390100e+06,
	    7.88027700e+06, 7.64807500e+06, 7.25954850e+06, 6.99782700e+06, 6.54252700e+06, 6.29859250e+06, 6.06362300e+06, 5.84142800e+06, 5.36775850e+06,
	    4.89408950e+06, 4.68576750e+06, 4.49476350e+06, 4.38545700e+06, 4.21975200e+06, 4.09634625e+06, 4.15417150e+06, 4.04259400e+06, 3.96940425e+06,
	    3.89933150e+06, 3.81049975e+06, 3.72820650e+06, 3.66312475e+06, 3.59110975e+06, 3.52811250e+06, 3.46507475e+06, 3.40187575e+06, 3.33257625e+06,
	    3.26777575e+06, 3.15637975e+06, 2.95494125e+06, 2.92330350e+06, 2.87439675e+06, 2.80228975e+06, 2.73674900e+06, 2.68266550e+06

	};

	// Interpolate to find the lifetime of a star
	auto const &x_arr = interp_mass_star;
	auto const &y_arr = interp_death_time;
	const double mass_in_Msun = mass_star / C::M_solar;
	AMREX_ASSERT(mass_in_Msun >= 0.0);
	amrex::Real death_time = interpolate_value<BoundaryPolicy::Clamp>(mass_in_Msun, x_arr.data(), y_arr.data(), AGE_ARR_SIZE) * YR_TO_SEC; // NOLINT
	return death_time;
}